---
title: "Analysis of models"
output: html_notebook
---
# Part 6: Analysis of models

# Load packages
```{r}
library(rethinking)
```

## Load data from previous modules
```{r}
# Load data from module 1
load("./data/d_sol_cw_w_in.RData")

# Load models and data from module 4
load("./data/model_out_final_real.RDATA")
load("./data/modeldata_out_final_real.RData")

```

## Run counterfactual simulation for solar radiation
```{r}
Sol_seq <-seq(from=0,to=5,length.out=1000)
sim_data <- data.frame(Sol_std = Sol_seq)
s1 <-sim(model_m1, data=sim_data, vars=c("T_S_std", "C_MW"))

# Inverse transform for plotting
Sol_actual <- Sol_seq * sd(d_sol_cw$Sol) + mean(d_sol_cw$Sol)

# Plot with de-standardised axes
windows(); plot(Sol_actual,colMeans(s1$C_MW),ylim=c(0,7),type="l", xlab=expression(Sol~"["*W/m^2*"]"), ylab=expression(C[MW] ~ "[" * pg/L * "]")) 
grid()
shade( apply(s1$C_MW,2,PI),Sol_actual)
mtext(expression("Regression of the effect of " * Sol * " on " * C[MW] * " using model " * m[1]))
```

```{r}
Sol_seq <-seq(from=0,to=5,length.out=1000)
sim_data <- data.frame(Sol_std = Sol_seq)
s4 <-sim(model_m4, data=sim_data, vars=c("T_S_std","C_MW", "W_std", "r_W_std"))

# Inverse transform for plotting
Sol_actual <- Sol_seq * sd(d_sol_cw$Sol) + mean(d_sol_cw$Sol)

# Plot with de-standardised axes
windows(); plot(Sol_actual,colMeans(s4$C_MW),ylim=c(0,7),type="l", xlab=expression(Sol~"["*W/m^2*"]"), ylab=expression("simulated " * C[MW] ~ "[" * pg/L * "]"))
grid()
shade( apply(s4$C_MW,2,PI),Sol_actual)
mtext(expression("Regression of the effect of " * Sol * " on " * C[MW] * " using model " * m[4]))

```


## Sample from posterior distributions
```{r}
dlist <- list(
    C_MW = d_sol_cw$C_MW,
    C_MW_w = d_sol_cw$C_MW,
    Sol_std = d_sol_cw$Sol_std,
    T_S_std = d_sol_cw$T_S_std
)

# Simple model
mu_po <- link(model_m1)
mu_po_mean <- apply(mu_po$mu, 2, mean)
mu_po_pi <- apply(mu_po$mu, 2, PI)

windows(); plot(mu_po_mean ~ dlist$C_MW, col=rangi2, ylim=range(mu_po_pi),
xlab = "Observed C_MW",ylab = "Predicted C_MW")
abline( a = 0, b = 1, lty = 2)
#for (i in 1:nrow(d_in))
#    lines(rep(dlist$C_MW[i],2), mu_po_pi[, i], col=rangi2)

# Partial model m2
mu_full_uinter <- link(model_m2)

mu_full_uinter_mean <- apply(mu_full_uinter$mu, 2, mean)
mu_full_uinter_pi <- apply(mu_full_uinter$mu, 2, PI)

windows(); plot(mu_full_uinter_mean ~ dlist$C_MW, col=rangi2, ylim=range(mu_full_uinter_pi),
xlab = "Observed C_MW",ylab = "Predicted C_MW")
abline( a = 0, b = 1, lty = 2)

# Full model with interaction
mu_full_inter <- link(model_m3)

mu_full_inter_mean <- apply(mu_full_inter$mu, 2, mean)
mu_full_inter_pi <- apply(mu_full_inter$mu, 2, PI)

windows(); plot(mu_full_inter_mean ~ dlist$C_MW, col=rangi2, ylim=range(mu_full_inter_pi),
xlab = "Observed C_MW",ylab = "Predicted C_MW")
abline( a = 0, b = 1, lty = 2)

# Full model with interaction and external confounders
mu_full_inter_wind <- link(model_m4)

mu_full_inter_wind_mean <- apply(mu_full_inter_wind$mu, 2, mean)
mu_full_inter_wind_pi <- apply(mu_full_inter_wind$mu, 2, PI)

windows(); plot(mu_full_inter_wind_mean ~ dlist$C_MW_w, col=rangi2, ylim=range(mu_full_inter_wind_pi),
xlab = "Observed C_MW",ylab = "Predicted C_MW")
abline( a = 0, b = 1, lty = 2)

postpredict <- list(
    m1_mean = mu_po_mean,
    m1_pi = mu_po_pi,
    m2_mean = mu_full_uinter_mean,
    m2_pi = mu_full_uinter_pi,
    m3_mean = mu_full_inter_mean,
    m3_pi = mu_full_inter_pi,
    m4_mean = mu_full_inter_wind_mean,
    m4_pi = mu_full_inter_wind_pi,
    C_MW = d_sol_cw$C_MW
)

save(postpredict, file = "./data/model_postpredictdist_out.RData")
```


## Calculate R-square values
```{r}
# Calculate R-square values
calculate_R_squared <- function(y_obs, model) {
  # Extract the posterior predictions from the model
  l <- link(model)
  y_mu <- apply(l$mu, 2, mean)
  
  # Calculate the total sum of squares (SS_tot)
  SS_tot <- sum((y_obs - mean(y_obs))^2)
  
  # Calculate the residual sum of squares (SS_res)
  SS_res <- sum((y_obs - y_mu)^2)
  
  # Calculate R-squared
  R_squared <- 1 - (SS_res / SS_tot)
  
  # Return the R-squared value
  return(R_squared)
}

y_obs <- d_sol_cw$C_MW

R_squared_m1 <- calculate_R_squared(y_obs, model_m1)
R_squared_m2 <- calculate_R_squared(y_obs, model_m2)
R_squared_m3 <- calculate_R_squared(y_obs, model_m3)
R_squared_m4 <- calculate_R_squared(y_obs, model_m4)

print(R_squared_m1)
print(R_squared_m2)
print(R_squared_m3)
print(R_squared_m4)
```